﻿/* Copyright 2007-2008 dnAnalytics Project.
 *
 * Contributors to this file:
 * Patrick van der Velde
 * Marcus Cuda
 *
 * This file is part of dnAnalytics.  dnAnalytics is licensed under the 
 * Microsoft Public License. See License.txt for a complete copy of the
 * license.
 */
using System;
using dnAnalytics.Properties;

namespace dnAnalytics.LinearAlgebra.Decomposition
{
    internal abstract class AbstractCholesky
    {
        private readonly Matrix mFactor;
        private readonly int mOrder;
        private bool mComputed;
        private bool mIsPositiveDefinite;
        private double mDeterminant = double.MinValue;

        protected AbstractCholesky(Matrix matrix)
        {
            mOrder = matrix.Rows;
            mFactor = matrix.Clone();
        }

        public bool IsPositiveDefinite()
        {
            Compute();
            return mIsPositiveDefinite;
        }

        public Matrix Factor()
        {
            if (!IsPositiveDefinite())
            {
                throw new NotPositiveDefiniteException();
            }

            return mFactor.Clone();
        }

        public double Determinant()
        {
            if (!IsPositiveDefinite())
            {
                throw new NotPositiveDefiniteException();
            }

            if (mDeterminant == double.MinValue)
            {
                lock (mFactor)
                {
                    mDeterminant = 1.0;
                    int rows = mOrder;
                    for (int j = 0; j < rows; j++)
                    {
                        double value = mFactor.ValueAt(j, j);
                        mDeterminant *= (value*value);
                    }
                }
            }
            return mDeterminant;
        }

        public Matrix Solve(Matrix input)
        {
            Matrix result = mFactor.CreateMatrix(mFactor.Columns, input.Columns);
            Solve(input, result);
            return result;
        }

        public void Solve(Matrix input, Matrix result)
        {
            if (input == null)
            {
                throw new ArgumentNullException("input");
            }
            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if (mFactor.Rows != input.Rows)
            {
                throw new NotConformableException("input", Resources.ParameterNotConformable);
            }

            if (mFactor.Columns != result.Rows)
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }

            if (input.Columns != result.Columns)
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }

            if (!IsPositiveDefinite())
            {
                throw new NotPositiveDefiniteException();
            }

            input.CopyTo(result);
            DoSolve(mFactor, result);
        }

        public Vector Solve(Vector input)
        {
            Vector result = mFactor.CreateVector(mFactor.Columns);
            Solve(input, result);
            return result;
        }

        public void Solve(Vector input, Vector result)
        {
            if (input == null)
            {
                throw new ArgumentNullException("input");
            }
            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if (mFactor.Rows != input.Count)
            {
                throw new NotConformableException("input", Resources.ParameterNotConformable);
            }

            if (mFactor.Columns != result.Count)
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }

            if (!IsPositiveDefinite())
            {
                throw new NotPositiveDefiniteException();
            }

            input.CopyTo(result);
            DoSolve(mFactor, result);
        }

   
        private void Compute()
        {
            if (!mComputed)
            {
                mIsPositiveDefinite = DoCompute(mFactor, mOrder);
                mComputed = true;
            }
        }

        protected abstract bool DoCompute(Matrix data, int order);

        protected abstract void DoSolve(Matrix factor, Matrix result);

        protected abstract void DoSolve(Matrix factor, Vector result);
    }
}
